/*
 *  component.h
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#ifndef _component_h_
#define _component_h_

#include <vector>
#include <string>
#include <map>
#include <algorithm>
#include <Eigen/Dense>

#include "point.h"
#include "MATERIALS/material.h"
#include "MATH/sparsematrixslu.h"
#include "commonutil.h"
#include "MATH/TENSOR/tensor.h"
#include "MATH/TENSOR/xtensor.h"
#include "MATH/TENSOR/tutils.h"


class Body;

using namespace std;
using namespace Eigen;

class Component {
public:
	int dim;					/** Dimension, from the parent body */
	int numDofs;				/** Number of dofs */
	int numDofsInertia;
	vector<string> pointNames;						/** Names of the points belonging to the component */
	map<string,Point*,less<string> > mPoints;		/** Pointers to points belonging to the component */
	// WARNING: pointNames and mPoints have different ordering
	vector<int> localGlob;		/** Local-global numbering */
	vector<int> localGlobInertia;		/** Local-global numbering for dofs contributing inertia */
	double time;				/** Time (local) */
	double deltaT;
	string matName;
	vector<double> q;			/** Local position dofs */
	vector<double> dq;			/** Local velocity dofs */
	vector<double> qC;			/** Local converged position */
	vector<double> dqC;			/** Local converged velocity */
	//vector<double> Fi;			/** Local contribution to the residual, inertia */
	VectorXd Fi;			/** Local contribution to the residual, inertia */
	VectorXd Fq;			/** Local contribution to the residual, function of q and dq */
	//vector<double> Ft;			/** Local contribution to the residual, function of t */
	VectorXd Ft;			/** Local contribution to the residual, function of t */
	//vector< vector<double> > Ki;		/** Local contribution to the tangent, inertia */
	MatrixXd Ki;		/** Local contribution to the tangent, inertia */
	MatrixXd Kq;		/** Local contribution to the tangent from Fq, derivative of q  */
	//vector< vector<double> > Kdq;		/** Local contribution to the tangent from Fq, derivative of dq */
	MatrixXd Kdq;		/** Local contribution to the tangent from Fq, derivative of dq */
	//vector< vector<double> > M;		/** Local contribution to the mass matrix */
	MatrixXd M;				/** Local contribution to the mass matrix */
	double perturbation;			/** For numerical differentiation */
	int integratorType;
	double thetab;		// Ambient temperature
	bitset<4> dofT;		/** Bitset for additive update dofs that belongs to the component */
	vector<int> globalTNum;	/** Global numbers of component dofs */
	vector<double> vLagMul;	// Lagrange multipliers
	bool augLagFlag, augConvergence;
	vector<double> Phi;		// Constraint vector
	vector<vector<double> > dPhi;	// Jacobian of constraint vector
	vector<double> PhiC;		// Coverged Constraint vector
	double penalty;
	vector<double> FqLag;
	vector<double> FqPenalty;
	double FqRatio;
	vector<double> MLumped;
	//vector< vector<double> > MSmooth;
	//vector<vector<double> > FSmooth;
	MatrixXd MSmooth;
	MatrixXd FSmooth;

public:
	Component();
	Component(vector<string>, int);
	Component(vector<string>, string, int);
	virtual ~Component(){}

	vector<string>& getPointNames();
	void setMPoints(map<string,Point*,less<string> >);
	void reportContent();
	void fillResidual(vector<double>& res_in_in,
					  vector<double>& q_in,
  					  vector<double>& dq_in,
					  const double& t, const int&, const double&, const int&, const double&);

	void fillResidual(	vector<double>& residue,
						const vector<double>& q,
						const vector<double>& qdot,
						const vector<double>& qddot,
						const double& t,
						const int&,
						const double&, const int&,
//JC						const int&, const vector<double>&, const vector<double>&,
  						const int&,
						const double&
					);

	void fillResFi(	vector<double>& residue,
				   const vector<double>& q,
				   const vector<double>& qddot,
				   const int&
				   );
	void fillResFq(	vector<double>& residue,
				   const vector<double>& q,
				   const vector<double>& qdot,
				   const int&,
   				   const int&,
				   const double&
				   );
	void fillResFt(	vector<double>& residue,
				   const double& time,
				   const int&,
				   const double&, const int&
				   );		
	
	void fillTangent(sparseMatrix&, vector<vector<double> >&, vector<vector<double> >&, vector<vector<double> >&,
					 vector<int>&, vector<int>&,
					 vector<double>&, vector<double>&,
					 const double& time, const int&);
	
	void fillTangent(sparseMatrix&, vector<vector<double> >&, vector<vector<double> >&, vector<vector<double> >&,
					 vector<int>&, vector<int>&,
					 const vector<double>& q,
					 const vector<double>& qdot,
					 const vector<double>& qddot,
					 double partial_qdot,
					 double partial_qddot,
					 const double& t,
					 const int&,
					 const int&,
					 const double&
					 );
	
	void fillMassMatrix(sparseMatrix& Ma, const vector<double>& q, const double& t);
	void fillMassMatrix(vector<double>& vMa, const vector<double>& q, const double& t);
	
	void fillForce(vector<double>& f,
						const vector<double>& q,
						const vector<double>& qdot,
						const double& t,
						const double&, const int&
					);

	void getLocalDofs(const vector<double>&);							/** Receive the global vector and extract local values */
	void getLocalDofs(const vector<double>&, const vector<double>&);	/** Receive the global vector and extract local values */
	void setLocalGlob();												/** Set local-global numbering */
	void setSparsePattern(vector< vector<int> >&);
	void setSparsePatternInertia(vector< vector<int> >&);
	void changePointName(string&, string&);
	void calcKqNumeric();
	void update(const int&, const vector<double>&, const vector<double>&);
	void setTemperature(double& th);
	int setTRNumbers(int&);
	void updateDofs(vector<double>&, vector<double>&);
	int setInertiaNumbers(int, vector<int>&);
	void updateLagMul(void);
	void checkConvLagMul(bool&);
	void calcProjMatrix(sparseMatrix&, double&);
	void calcCon(vector<double> &, vector<double>&, vector<double>&, vector<double>&, double&);
	void fillPhiq(vector<vector<double> >&, int&);
	void setPointNames(vector<string>&);
	void fillMSmooth(sparseMatrix& Ma, const vector<double>& q, const double& time );
	void fillFSmooth(vector<vector<double> >&, const vector<double>&, const vector<double>&, const double&);

	int checkComponents(int);
	
	
	/** Virtual functions */
	virtual void setDofs(void)=0;
	virtual void initializeGlobalDofs(vector<double> &, vector<double> &, const double &)=0;			/** Set initial dof values */
	virtual void calcFi(const vector<double>&)=0;
	virtual void calcFq(const int&)=0;
	virtual void calcFt(const double&, const double&, const int&)=0;
	virtual void calcKi(double&)=0;
	virtual void calcKq(int)=0;
	virtual void calcKdq(double&)=0;
	virtual void calcM(const vector<double>&, const double&)=0;
	virtual void calcMLumped(const vector<double>&, const double&)=0;
	virtual void setMatAtComp(Material*)=0;
	virtual string getMatName()=0;			/** Get the material name */
	virtual void initializeCompRes(const bool&, const double&, const int&, const bool)=0;
//	virtual void initializeCompTang(const bool&, const bool&)=0;
	virtual void initializeCompTang(const bool&)=0;
	virtual void calcEner(double&, double&, double&, double&, double&, double&, double&, int&, vector<double>&, vector<double>&)=0;
	virtual void calcEntropy(double&, double&)=0;
	virtual void calcAM(vector<double>&, vector<double>&, vector<double>&)=0;
	virtual void updateComp(void)=0;
	virtual void setInitialConditions(vector<double>&, vector<double>&)=0;
	virtual double getTemperature(void)=0;
	virtual void setIniDofs(vector<double>&)=0;
	virtual void calcPhi(const vector<double>&)=0;
	virtual void calcDPhi(const vector<double>&)=0;
	virtual void calcMSmooth(const vector<double>&, const double&)=0;
	virtual void calcFSmooth(const vector<double>&, const vector<double>&, const double&)=0;
	virtual int checkState(void)=0;
	
	/** Friend functions */
	friend void printSTLVector(vector<int>& v);
	friend void printSTLVector(vector<double>& v);
	friend void printSTLMatrix(vector< vector<int> >& v);
	friend void printSTLMatrix(vector< vector<double> >& v);
	friend double lTwoNorm(const vector<double>& v1);
	friend double lTwoNormSq(const vector<double>& v1);

};
#endif
